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ABSTRACT 


Finite-difference  methods  have  come  into  wide  use  for  solving 
special  problems  Including  transient-heat  conduction.  Dusinberre'1' 
has  ably  presented  the  possibilities  of  finite-difference  methods. 

The  success  of  most  such  methods  depends  on  the  existence  of  a 
certain  degree  of  uniformity  of  behavior  of  the  temperature  over  the 
finite  intervals  of  both  space  and  time  selected  for  the  computation 
process.  In  some  cases,  however,  this  required  uniformity  constitutes 
a  handicap  since  temperatures  are  changing  so  rapidly  that  Incon¬ 
veniently  short  time  intervals  have  to  be  chosen.  This  paper  repre¬ 
sents  an  effort  to  develop  a  finite-difference  method  free  from  the 
foregoing  defect. 


"Numerical  Analysis  of  Heat  Flow,"  by  C-. 
McQraw-Hill  Book  Company,  Inc.,  New  York.  N.  Y,, 


M,  Dusinberre, 
1949. 
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base  of  natural  logarithms,  2.718 
dummy  variables 

thermal  dlffuslvity,  ft^/hr 

dummy  time  variable,  hr 

function  describing  the  ambient  temperature,  °F 

Numerical  subscripts  apolied  to  temperatures,  as  in 
Ti»  To,  refer  to  temperatures  located  at  x  =  Ax, 

x  =  0,  and  x  =  "<Ax,  etc.,  at  the  time  t  =  0. 

Positive  and  negative  superscripts  applied  to  temperatures, 
mean  that  the  temperatures  are  to  be  evaluated  at  the  plus  and 
minus  sides  of  the  points  in  question. 

Introduction; 

In  recent  years  finite-difference  methods  have  been  used 
to  solve  a  vast  number  of  special  problems  involving  transient 
heat  conduction.  In  flexibility  of  use  and  simplicity  of  con¬ 
cept  these  methods  excel  those  of  classical  mathematical  analysis, 
and.  Indeed,  on  many  occasions  they  are  not  to  be  regarded  as 
substitutes  for  more  precise  methods,  but  as  the  only  possible 
methods  to  use.  The  possibilities  of  finite-difference  methods 
for  many  problems  akin  to  those  considered  in  this  paper  are 
well  presented  in  the  book  by  Dusinberre1. 

The  success  of  most  finite-difference  methods  depends  on 
the  existence  of  a  certain  degree  of  uniformity  of  behavior  of 
the  temperature  over  the  finite  intervals  of  both  space  and  time 
selected  for  the  computation  process.  There  are,  however,  a 
number  of  occasions  when  this  requirement  of  uniformity  is  a 
handicap,  since  temperatures  are  changing  so  rapidly  that  incon- 


veniently-short  tins  intervals  have  to  be  chosen.  This  awkward 
feature  usually  arises  near  the  boundary  of  the  computation 
region;  for  example,  it  may  arise  near  the  surface  of  a  casting 
during  quenching. 

The  present  investigation  represents  an  effort  to  develop 
a  finite-difference  method  free  from  the  foregoing  defect. 
Formulas  are  derived  which  do  not  imply  a  uniformity  of  behavior 
with  respect  to  time.  Within  the  interior  of  a  solid,  these 
formulas  reduce,  in  general,  to  those  obtainable  by  the  simpler 
technique  of  heat  balances.  However,  near  a  convective  heat- 
transfer  surface  they  do  not  reduce  to  earlier  formulas.  In 
this  region  they  possess  greater  potentiality,  in  that  they 
will  handle  with  uniform  precision  cases  of  variable,  and  even 
discontinuous,  ambient  temperature,  with  the  heat-transfer 
coefficient  ranging  from  zero  (insulation)  to  infinity  (perfect 
contact . ) 

In  actual  practice,  the  new  formulas  merely  introduce  new 
weighting  factors  into  the  standard  finite-difference  equations. 
A  numerical  table  of  such  weighting  factors  is  given  for  the 
temperatures  on,  and  adjacent  to,  a  convective  heat-transfer 
surface  when  the  space  and  time  intervals  are  chosen  to  conform 
with  the  Binder-Schmldt  selection  of  M  =  (Ax)^/kA  t  =  2- 
Derivation  of  Formulas  for  the  Infinite  Medium; 

Formulas  will  here  be  derived  which  are  appropriate  for 
use  in  computing  the  transient  conduction  of  heat  in  an  infin¬ 
ite  medium  of  uniform,  constant  properties.  The  differential 
equation  to  be  applied  is  as  follows: 

<LT-  K 

at  “ 

This  is  a  linear  differential  equation,  and  the  response  at 


(1) 
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some  time  "t"  la  linearly  related  to  the  temperature  dlstribu- 

tion  Input  "f(x)"  at  time  zero  by  the  solution  given  below  . 

/  r°°  , 

T(x,t)  •JySZr  Jf(x)  e  d*  (2) 

~Oo 

At  the  point  x=0  this  last  equation  reduces  to: 

-JL2 


T(0‘ *>  =  £^Tt 


Now  let  the  contribution  to  T(0,t)  originating  within  the 
region  x>0  be  denoted  by  C-^O.t).  This  contribution  can  be 
written  in  the  following  manners:  _ 

OO  y/^ 

i  /  ^ 


W>-fi A' 


'  '  '  zyirxt  J - 

0r:  °co 

C.(o,t)  =  +  (j(ztVKt)t  clir  (5) 

O 

The  integral  in  eq.  5  can  be  successively  integrated  by  parts 
to  give: 

C  (o,  t)  =  jf^f<(o)(Zt/Xt)i'erMo)  +  7 

'  **  L  r=  o  oo  rt4i  „u  ~] 

+  ffc'rtt)  /  (2f  VmJ)  J* J  (6) 

where  the  lnerfc( ^ )  are  the  iterated  error  functions  defined 


i'erfc  tk)  —  Yrr  6 


and: 


/  "erfi c  ft)  - 


The  presumption  in  writing  eq.  6  is  that  the  first  "n" 
derivatives  of  f(x*)  are  continuous  on  the  plus  side  of  x'=0 
These  derivatives  may,  or  may  not,  be  continuous  through 
x'=  0.  If  they  are  continuous,  then  when  the  contribution 


"Conduction  of  Heat  In  Solids,'  by  Carslaw  and  Jaeger1 
Oxford  University  Press,  19^7 • 


GgCOjt)  is  added  to  C^(0,t),  all  odd  temperature  derivatives 
cancel,  leaving  ths  following  result  for  T(C,t). 


r-o  2  o  (9) 

When  the  last  integral  in  ec.  9,  which  is  the  error  term, 


is  neglected,  and  when  use  is  made  of  the  fact  that: 


i%rh  (P) 


eq.  9  then  reduces  to 


T(c}t)~ 


n>(Kt)rf”(o) 


(10) 


(ID 


Thus  when  the  temperature  distribution  in  an  infinite  medium 
conforms  at  t=G  with  some  polynomial  in  '  x"  ,  eq.  11  gives  the 
exact  answer  for  all  subsequent  time.  In  terms  of  the  modulus 


"M"  this  equation  becomes: 

n 


2r  ,2r. 


where : 


T(0,  t)  =  2:  M  [(Ax)  *-~0)] 


(12) 


t-o  =  At 


(13) 


Equation  13  is  still  not  a  convenient  expression  to  use 
numerically.  An  Improvement  can  be  made  by  evaluating  the 
derivatives  in  this  expression  in  terms  of  discretely-spaced 
temperatures  with  the  uniform  Interval,  /^x.  Thus,  for  a 


second -degree  polynomial: 

(Axf  f°(o)  =  Tc  (1A) 

(Axf  fz(o)  =  T.-ZT.+T, 

For  higher-degree  polynomials,  the  second  and  higher  derivatives 


take  on  more  complicated,  but  similar,  form.  Substitution  of 
eqs.  14  into  eq.  12  gives  the  standard  finite-difference  form, 


widely  used  in  heat-conduction  studies.  Thus: 

T(o,Ai)  =  T.  +Ji(Tl-ZT0+Tl) 


(15) 


Case  of  M  =  rr : 


Equation  15  can  be  deduced  very  much  more  quickly  by 


direct  use  of  heat  balances.  In  this  case  the  ability  of  the 
present,  more  elaborate  analysis  to  accomodate  non-uniform 
time  behavior  is  not  made  evident.  To  bring  this  ability  into 
evidence,  let  it  be  supposed  now  that  at  time  zero  neither 
the  temperature  nor  its  derivatives  is  continuous  through 
x=0.  Such  a  situation  can  arise  physically  when  two  plates 
of  similar  material,  but  dissimilar  temperatures,  are  suddenly 
brought  into  good  thermal  contact.  It  obviously 
includes  as  a  speolal  case  the  more  usual  situation  analyzed 
above. 

In  the  present,  more  general  case,  0^(0, t)  and  C2(0,t) 
must  be  evaluated  separately.  Now  let  it  be  assumed  that  the 
temperature  distribution  at  t=0  for  x>0  can  be  well  represented 
by  a  second -degree  polynomial.  Furthermore,  let  the  tempera¬ 
tures  on  the  plus  side  of  varloue  stations  be  identified  with 
"plus'*  superscripts,  and  the  temperatures  on  the  minus  side 
with  "minus*  superscripts.  Then: 

t°(+a)=  Tc  + 

(&x)f,(+0)=-k(+T+-3T,-X)  (i6) 

(Ax)Y(+0)=  Tz-ZX  +  V 

Since  ierfc(O)  =  1//tt  ,  use  of  eq.  16  in  eq.  6  gives*. 

C,(o,t)=  [z  (tyfif+M')'*'  ~Vf)  *  ^2  {~M  ~vWrr)J  (17) 

This  last  expression  is  valid  regardless  of  the  temperature 
distribution  for  x<0;  that  is,  regardless  of  the  magnitude  of 
the  time  derivatives  induced  by  discontinuities  at  x=0. 

Inspection  of  eq.  17  shows  that  a  rather  remarkable 
simplification  occurs  if  M  -  n.  In  this  event,  only  the  tem¬ 
peratures  Tq  and  T-^  need  to  be  known  in  order  to  obtain  exact 


results  for  an  Initial  quadratic  temperature  profile.  If 
the  profile  for  x<0  Is  also  quadratic,  though  different, 
addition  of  the  contributions  C-^  and  leads  again  to  eq.  15, 
provided  the  temperature  Tq  in  that  equation  be  Interpreted  as: 

X~i(V+X)  (18) 

Thus  the  new  method  of  deducing  the  finite-difference  equations 
has  brought  out  the  unique  property  of  M  -  tt; 


namely,  that  it  can  accomodate  space  discontinuities  in  the 
temperature  and  its  derivatives,  if  these  discontinuities 


occur  at  a  central  grid  point. 

Although  the  interpretation  of  Tq  according  to  eq.  18 
makes  eq.  15  for  M  =  tt  highly  accurate  when  a  temperature- 
jump  occurs  at  the  central  grid  point,  there  remains  the 
problem  of  how  to  weight  the  temperatures  at  station  1,  say, 

If  a  temperature  Jump  occurs  there,  instead.  As  before,  let 
the  temperature  profiles  to  the  right  and  left  of  station  1 
be  quadratic,  though  different.  Then  C2(0,t)  can  immediately 
be  written  as: 

Q  a?) 

On  the  other  hand,  C^(0,t)  requires  special  treatment. 

Let  f  (x1 )  be  the  smooth,  or  analytic,  continuation 
Into  the  region  x'>Ax  of  the  actual  temperature  profile 


existing  in  the  region  x'<Ax.  Then  0^0,0  can  be  written 

c, (°> *)  +  {!%') e** dx'+j 


+  fyTTt  I  {!(*')- fbc)}  6+**  dx 

Ax 

The  first  two  integrals  in  this  last  equation  can  be  summed 

iR-iZ  +  iTT] 


(20) 


for  M  =  tr  to  give: 


The  third  integral  might  be  integrated  by  parts,  as  on  earlier 

occasions,  to  yield  a  series.  However,  since  such  a  procedure 

would  complicate  any  formula  by  introducing  temperatures 

beyond  T^,  it  is  not  followed.  Instead,  the  functional 

■if" 

difference  f (x1 )  -  f  (x* )  is  treated  as  stationary  compared 
with  the  fast -attenuating  exponential,  and  the  last  integral 
is  approximated  by  the  following  expression: 


( T*-T~)i  (£)  =  o.  1053  (if-  7 f) 

Combination  of  the  foregoing  results  gives  the  complete 
expression  for  C^;  i.  e., 

*  #  fj  +  ojossft- 1 fj 

Since  0.1053  n  =  0.331  =  1/3,  the  sum  of  C^+Cg  can  be 

written  with  high  accuracy  as:  2T~+  T+ 

,=  1.  TL-2Z+-Lr' 


(21) 


(22) 


T(o,m)=  4  +  a.-* _ f.  (23) 

Comparison  of  eq.  23  with  eq.  15  shows  that  the  standard  form 
of  eq.  15  can  be  retained  if,  when  predicting  T(0,At)  with  a 
temperature  jump  at  station  1  (x  =Ax),  the  "inside"  temperature 
at  the  jump  is  weighted  twice  as  heavily  as  the  "outside”  tem¬ 


perature.  Or,  in  other  words,  in  eq.  15  is  to  be  interpreted 
as:  Tx  =  2  T1  *  T1 


(24a) 


Put  alternatively,  the  following  extended  form  of  eq.  15  yields 
excellent  accuracy  when  M  =  rr  and  when  the  temperature  profile 
can  be  represented  by  second -degree  curves  in  the  intervals: 
-Ax  ;  -Ax^x  <  0  J  0<rx  <  Ax  ;  X  >■  Ax 

Thus!  T(o,tA)-±(£*  X) 

To  illustrate  the  practical  use  of  the  foregoing  rules, 
let  us  consider  the  following  example.  "A  semi-infinite  slab 


(x>0)  at  a  temperature  of  100  degrees  is  suddenly  brought  into 
perfect  thermal  contact  with  a  second  semi-inf Inlte  slab  (x<0) 
at  a  f  amperature  of  -100  degrees.  The  temperature  distribution 
for  all  time  for  x>C  Is  desirod." 

The  above  example  13  solved  in  ref.  1  (p.  121)  for 
M  =  2,  3  and  A.  Table  I  of  this  paper  shows  the  numerical  results 

obtained  when  M  =  n.  For  the  chosen  mathematical  model, 
exact  values  are  shown  in  parentheses.  At  no  tabulated 
point  does  the  absolute  error  of  the  finite-difference  method 
exceed  3«2&. 

TABLE  I 


-  Ax 

0 

Ax 

2  Ax 

3  Ax 

b  Ax 

0 

-100 

IOC 

100 

100 

100 

At 

0 

78.7 

(79.0) 

100 

(98.8) 

100 

100 

2  At 

0 

60.5 

(62.4) 

93.3 

(92.3) 

100 

100 

3  At 

0 

51.7 

(53.0) 

CD  00 

•  * 

ro  O 

97.9 

(97.0) 

100 

It  is 

interesting 

to 

note  that 

the  temperature 

at  x=0 

la,  for  the  purposes  of  computing  future  temperatures  at  the 
same  point,  taken  as  zero  right  from  the  start.  However, 
for  purposes  of  computing  future  temperatures  at  station  1, 
the  Initial  temperature  at  x=0  Is,  by  eq.  24,  to  be  taken 
as:  — — °°y-  (~10Q)  =  33.3  deg. 

The  major  source  of  improvement  of  the  present  computational 
accuracy  over  the  cited  examples  In  ref.  1  lies  in  the  treat¬ 
ment  of  temperatures  at  a  point  of  discontinuity. 

If  automatic  computing  machinery  Is  to  be  used  for  the 
finite-difference  computations,  the  selection  of  M  =  tt  should 
introduce  negligible  inconvenience.  However,  for  hand  compu- 
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fcation  the  use  of  M  =  3,  a  value  quite  close  to  rr,  would 
appear  to  be  preferable  because  of  the  simpler  arithmetical 
manipulations  required.  The  rules  given  by  eqs.  18  and  24 
should  be  retained.  In  the  foregoing  example,  use  of  M  =  3 
with  these  rules  increases  the  maximum  error  by  only  C.3%. 

Formulas  for  the  Neighborhood  of  a  Convective  Heat-Transfer 

Surface: 

The  formulas  obtained  in  the  previous  section  are  valid 
for  the  infinite  medium,  or  for  finite  regions  which  can  be 
mimicked  by  an  infinite  medium  through  the  use  of  superposition 
of  symmetric  and  anti-symmetric  temperature  distributions.  In 
the  general  case  of  heat  convection  from  a  surface,  however, 
the  heat  transfer  coefficient  is  usually  neither  so  small  that 
heat  transfer  can  be  neglected,  nor  so  large  that  perfect 
thermal  contact  can  be  assumed.  This  general  case  does  not, 
unfortunately,  lend  itself  readily  to  the  superposition  technique, 
and  special  formulas  are  required.  Suitable  formulas  of  high 
accuracy  will  be  given  in  this  section.  Their  detailed  derivation 
is  given  in  the  Appendix,  and,  because  no  new  principles  are 
involved,  it  will  suffice  here  to  summarize  and  illustrate  the 
results. 

The  short-term  behavior  of  all  slabs  of  finite  thickness 
is,  with  respect  to  changes  at  their  surfaces,  like  that  of 
corresponding  semi-infinite  slabs.  Accordingly,  the  results 
obtained  for  the  semi-infinite  slab  whose  surface  is  exposed 
to  convective  heat-transf er,  can  also  be  used  for  the  finite 
slab,  so  long  as  At  for  the  time  interval  of  computation  is 
not  too  large.  Consider,  therefore,  a  semi-infinite  solid 
medium  having  uniform,  constant  properties.  Within  the  solid 
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the  temperatures  must  obey  Fourier's  equation  (eq.  1).  At 
the  surface  of  the  solid,  cooling  takes  place  according  to 
Newton's  "Law"  of  Cooling;  1.  e.. 


(25) 


Now  let  the  temperature  within  the  solid  at  t  =  0  be 
expressible  as  a  second -degree  polynomial  in  "x"  for  all 
"x"  beyond  the  surface.  Further,  let  the  ambient  temperature 
between  t  =  +0  and  t  =  At  -  0  be  a  linear  function  of  time. 
Then  because  of  the  linearity  of  the  governing  differential 
equation  and  its  boundary  conditions,  the  temperatures  T(C,At) 
and  T(Ax,At)  can  be  linearly  expressed  in  terms  of  T*,  T^, 
T2,  Ta(+0)  and  Tg(At  -  0).  The  weighting  factors  for  the 
various  temperatures  are  arrived  at  in  a  manner  similar  to 
that  used  for  the  case  of  the  infinite  medium.  The  results 
are:  T(0,  At)  =  AqtJ  +  +  C0T2 

D0Ta(+0)  ♦  E^Ta(At  -  0)  (26) 

T(  Ax,  At)  s  AxtJ  +  +  OjTg 

DlTa(+0)  +  ExTa(At  -  0)  (27) 

The  coefficients  defined  by  eqs.  26  and  27  are  given  by  the 
following  formulas: 


—  -N  F*"  ~  ^  p*  -  m  p*  4.  1 

J  NF1  2  a  N  *3  *  M 

J  =  1,2 

(28) 

j  =  2FV  K  + 

«l 

(29) 

PT  1 

J  “  21  +  M  ”  N  F3 

If 

(30) 

j  =  *  NM  h 

II 

(3D 

,  =  NM  F, 

J  3 

tl 

(32) 

In  these  equations  "N"  is  a  Nusselt  number  defined  by: 

N  =  h(Ax)/k 


(33) 


1 


The  various  F's  are  dependent  on  the  choice  of  "J" , 
to  be  calculated  by  the  following  relations: 

-JtM  ' 


where: 


and  are 

(34) 

(35) 


fr\^[e,fc(^)-Fc]  (36> 

'36) 

The  coefficients  in  eq.  2 6  and  27  are  functions  of 
M,  N  and  the  position  parameter,  j.  Their  calculation  in¬ 
volves  a  fair  amount  of  numerical  work.  However,  for  specific 
and  widely-used  values  of  M, tables  of  these  coefficients  can 
be  prepared  for  universal  use.  One  such  table,  for  M  =  2,  is 
given  in  the  next  section.  (Table  II)  It  occupies  little 
space,  yet  is  suitable  for  linear  interpolation  over  the 
entire  range  of  possible  values  of  the  heat  transfer  coefficient. 
When  such  a  table  is  available,  use  of  the  new  coefficients 
is  very  straightforward.  Questions  of  stability  do  not  arise, 
and  discontinuities  of  temperature  in  both  space  and  time  at 
the  surface  are  handled  automatically. 

F  t  grid  points  more  than  distance  Ax  from  the  surface,  the 
standard  finite-difference  formula  appropriate  to  the  chosen  M 
should  be  used.  (See  eq.  15).  At  a  sacrifice  of  accuracy, 
this  formula  can  also  be  used  to  compute  the  temperature  history 
at  x  =  ^x. 
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The  Special  Case  of  M  =  2: 

To  illustrate  the  capabilities  of  the  new  finite-difference 
formulas,  a  table  for  M  =  2  will  now  be  given,  and  used  to 
solve  typical  problems.  Table  II  gives  the  needed  coefficients. 
Equations  26  and  27  are  repeated,  making  use  of  the  table  nearly 
self-explanatory.  The  argument  to  be  used  in  entering  the 
table  is  l/(N+l),  which  is  +ha  rntio  of  the  surface  resistance 

to  the  sum  of  the  surface  res*,  stares  rlus  the  resistance  of 
a  slab  of  thickness ,  Ax.  A  simple  check  which  all  such 
tables  must  satisfy  is  th s>;  the  of  'he  coefficients  for 
any  g!  ran  rust  '  unity.  Th't  -,tis  statement  must  be  true 
can  be  seen  from  the  fact  that  if  all  temperatures  within  the 
solid  are  the  same  as  the  constant  ancient  temperature,  the 
temperatures  at  x  =  0  and  x  =  Ax  must  be  the  same  at  the  end 
of  time  At  as  at  the  beginning.  When  the  ambient  temperature 
is  constant,  very  often  it  can  be  used  as  the  datum  temperature 
(1.  e.f  taken  as  zero),  and  in  this  event  the  coefficlenta 
Dj  and  Ej  do  not  enter  the  computation. 

Consider  the  following  problem  to  illustrate  the  use  of 
the  foregoing  table.  "A  semi-infinite  medium  of  uniform,  con¬ 
stant  properties  is  everywhere  at  a  temperature  of  1000  degrees 
at  time  zero.  At  this  time  convective  cooling  commences  at 
its  exposed  surface  to  an  ambient  temperature  of  0  degrees.  The 
thermal  properties  of  the  medium  are  known,  as  well  as  the 
value  of  the  surface  heat-transfer  coefficient.  Find  the  tem¬ 
perature  history  within  the  slab." 

This  problem  is  solved  by  the  new  technique  by  selecting 
first  a  size  of  space  Interval  suitable  for  sampling  the  tsm- 


TABLE  II 


BOUNDARY  INFLUENCE  COEFFICIENTS  FCh  M  =  2 


T(0,4t)  =  A  T*  +  B  T,  +  C  T.  +  D  T  (+C)  +  ST  ( At  -  0) 
oo  oi  o  d  oa  oa 


I 

A 

0 

3o 

Co 

Do 

Eo 

N  +  1 

0 

© 

0 

0 

0 

1 

0.1 

0.0129 

0.0480 

0  .0267 

0.0672 

0.8452 

0.2 

0.0334 

0.1C80 

0.0474 

0.1087 

0 .70  25 

0.3 

C.0609 

0.1749 

0.0628 

0.1276 

C.5738 

0.4 

0.0935 

0.2438 

0.0743 

C . 1292 

0.4592 

0.5 

C . 1290 

0.3116 

0.0826 

0.1189 

0.3579 

0.6 

0.1654 

0.3765 

0.0888 

0.1C09 

0 . 2684 

0.7 

0.2017 

0.4375 

C  .0933 

0.0783 

C.1892 

0.8 

0 . 2370 

0.4943 

0.0967 

0.C531 

0.1189 

0.9 

0.2709 

0.5471 

0.0991 

0.0266 

C.C563 

1.0 

0.3032 

0.5957 

C.1C11 

0 

C 

1 

N  +  1 

A1 

B1 

C1 

D1 

E1 

0 

0 . 1074 

0 . 1507 

0.4246 

0 . 1 666 

0.1507 

0.1 

0.1256 

0.1802 

0.4248 

0.1491 

0.1203 

0.2 

0.1443 

0.2074 

0.4243 

0.1286 

0.0954 

0.3 

0.1624 

0.2315 

0.4234 

0.1077 

0.0750 

0.4 

0.1792 

0.2527 

0 . 4222 

0.0877 

0.0582 

0.5 

0 . 1944 

0.2710 

0.4212 

0.0692 

0  .0442 

0.6 

0.2081 

0.2869 

0.4201 

0.0524 

0.0325 

0.7 

0.2204 

0.3008 

0.4191 

0.0372 

0.0225 

0.8 

0.2314 

0.3130 

0.4182 

0.0235 

0.0139 

C.9 

0.2414 

0.3234 

0.4176 

C  .0114 

0.0062 

1.0 

0.2500 

C .3333 

0.4167 

0 

0 

perature  distribution  in  the  regions  of  interest.  The  time 
Interval  mus  t'r.n  be  ch;> -sen  sc  that  M  =  2.  Also,  from  the 
space-interval  selection,  the  surface  Nuaselt  Number,  N,  can 

be  calculated.  This  last  parameter  determines  the  coefficients 
which  are  read  from  Table  II.  In  the  present  case,  suppose 
N  =  1/2.  Table  III  gives  computed  results  for  six  time  Intervals 
Certain  exact  results  are  given  in  parentheses  to  provide  a 
gauge  of  the  computational  accuracy.  At  no  point  of  comparison 
does  the  error  of  the  finite-difference  process  exceed  0.1%. 

TABLE  III 


Ta 

0 

Ax 

2  Ax 

3  Ax 

4  Ax 

5  Ax 

6  Ax 

0 

0 

1000 

1000 

1000 

1000 

1000 

10CC 

1000 

At 

c 

699.2 

(699.2) 

932.4 

1C0C 

1000 

1000 

1000 

1000 

2  At 

0 

613.8 

(615.7) 

847.1 

966.2 

1000 

100C 

1000 

1000 

3At 

c 

558.9 

(562.6) 

789.2 

923.6 

983.1 

1000 

1000 

1000 

4  At 

0 

520.4 

(523.2) 

742.3 

886.2 

961.8 

991.6 

1000 

1000 

5  At 

0 

490.1 

(492.6) 

704.4 

852.0 

958.9 

98C.  9 

995.8 

1000 

6At 

0 

465.3 

(467.2) 

672.2  821.7 
(674.7) (822.9) 

916.4 

967.4 

990.5 

997-9 

To  assess  the  worth  of  the  present  adaptation  of  the 
method  of  finite-differences,  one  must  compare  it  with  alterna¬ 
tives.  For  example,  as  shown  in  ref.  1,  p.  129,  a  heat-balance 
at  the  surface,  made  on  the  assumption  of  constant  temperature 
gradients  throughout  one  time  Interval,  gives  coefficients 
equivalent  in  application  to  the  present  Aq ,  8^  and  DQ .  The 
formula  is  as  follows: 

T(0,At)  =  ^  Ta  4  M  T1  +  £  "  M(N+1i/  TC  (39) 

It  is  to  be  used  in  conjunction  with  eq.  15  for  all  interior 


points.  Dusinberre1  shows  that  stability  In  the  numerical  cal¬ 
culations  requires  that  M  be  greater  than  2(N+1)  In  eq.  39,  and 
greater  than  two  In  eq.  15-  (Thus  M-2  used  In  the  above  example 
Ip  at  the  limit  of  stability  of  eq.  15  and  beyond  the  limit  of 
stability  of  eq.  39.)  When  N  =  1/2,  M  =  4  meets  the  foregoing 
stability  criteria,  and  this  value  of  the  modulus  was  used 
with  eqs.  15  and  39  to  solve  the  example  problem.  Retention 
of  the  same  space  interval  meant  the  use  of  twice  as  many  time 
intervals  to  achieve  the  same  real  time;  that  is,  twice  as 
many  computation  points  were  required.  The  error  in  this  altern¬ 
ative  calculation  was  almost  uniformly  twice  as  great  as  in  the 
calculation  tabulated  in  Table  III. 

To  illustrate  further  the  use  of  the  table  of  coefficients, 
two  other  problems  will  be  solved  for  the  first  few  time  inter¬ 
vals.  As  a  first  illustration,  suppose  that  the  heat-transfer 
coefficient  in  the  problem  Just  solved  were  essentially  infinite 
Then  the  solution  of  the  problem  would  start  as  shown  in 
Table  IV. 


TL 

+0 

Table  IV 

Ax 

2Ax 

3  Ax 

"r  ~  1-1 

0 

0 

1000 

1000 

1000 

1000 

At 

0 

0 

683 

loco 

1000 

ZAt 

0 

0 

528 

842 

1000 

Table  V 


z 

ft 

+0 

A*. 

3  Ax. 

0 

u 

100 

200 

300 

400 

At 

0  *4 
** 

180 

217 

3co 

400 

2At 

214 

242 

308 

400 

As  a  second  example,  consider  a  ■ emi -Infinite  slab  having 
uniform  thermal  properties.  Let  the  exposed  surface  be  Insulated 
and  let  the  temperature  vary  linearly  with  distance  from  the 
exposed  surface.  The  calculations  begin  as  in  Table  V. 
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APPENDIX 

In  this  Appendix  detailed  derivations  are  given  for  the 
functions  and  formulas  useful  in  calculating  heat  conduction 
in  a  solid  near  a  convective  heat-transfer  surface. 

PROPERTIES  OF  THE  FUNCTIONS  Fn(x,  t,4)  and  Gn(x) 

The  Iterated  Error  Functions : 

The  Iterated  Error  Functions  are  defined  by  eqs.  7  and  8. 
They  are  tabulated  in  refs.  2.  From  eq.  7  differentiation 

gives: 

ferfc  (x)  —  **  erfc  (x)  (4o) 

dx 

This  last  relation  oan  be  used  to  give  meaning  to  iterated 

functions  with  Indices  less  than  minus  one  (-1).  Thus: 

i  *erfc  (x)  =  £  (41) 

o 

The  recursion  equation  satisfied  by  these  functions  is  : 

Zn  Cerfc  Oc)**  i*~zerfc(x)  -  Zx  i*  \ erfc(x )  (42) 

Equations  7  and  8  can  be  used  to  establish  the  validity  of 
eq.  42  for  all  11  n". 

4  2 

Definition  of  Fn(x,t,*^): 

The  functions  "Fn"  are  defined  by  the  recursion  equation: 


(43) 


(44) 


[(Z's/Kt)  ~  Fn] 

"1th:  S  4£* 

F0  *  b+*  erfc  {fki  +  Ai/Kt}  (44) 

Since  both  FQ  and  all  the  (2V*t)  *<-  "e  rfc  (^-t)  satisfy  the 
heat-conduction  equation,  it  follows  that: 

0  -  i  ft”  <45) 

Now  suppose  that,  for  a  moment: 

0*>  -  -  f„  («> 

then:  ,  r  *-i  /  v  .  “7 

1  etfc(FVKt^  “  Fx  J  (47) 

But  also :  - Fn  =£'[-(? FKt)  +  Fn_,  J  («> 

Therefore: 

7j?  =  ~  Fn -/  (49) 

Thus  eq.  49  Is  true  if  eq.  46  is  true.  But  likewise,  if 
eq.  49  is  true,  eq.  46  can  be  proved  to  be  true.  Finally, 
eq.  46  can  be  directly  verified  for  n  =  0.  Hence,  by  Induction, 
eq.  46  is  true  in  general. 

Also,  through  the  use  of  eqs.  45  and  46,  a  second  useful 
relation  can  be  found.  Thus: 

=  F  =  -L  (50) 


(45) 


(46) 


(47) 


(48) 


(49) 


p  -  J- 

r»-i  K 


'/»-/  "  K  Jt 

Recursion  Formula  for  the  G-  Functions: 


for  all  "n'1 


(51) 


A  second  type  of  function  appears  in  the  heat-conduction 


formulas  to  be  developed.  It  is  defined  by: 

Q-h  (x)  =  in  erf c  (x)  +  inerfc(-x) 

A  complementary  set  of  functions  is  defined  by: 

Hh  (x)  s  Perfc  (x)  -  i  "erfc  C-x) 


(52) 


(53) 
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The  recursion  equation  for  the  Iterated  Error  Functions  easily 

gives:  -  Cc  =  -Zx.  H„_,  (54) 


and: 


2n  Hh  -  Hn.z  CTh_t 


(55) 


When  values  of  "Hn"  from  eq.  54  are  substituted  into  eq.  55, 
the  following  recursion  equation  is  obtained  for  the 
Thus :  ,  a 

Q  =  - r-",-"-',  '  ~'"~J  (56) 


n 


'n+i  2n(n+i) 

The  first  few  functions  are  given  below. 

Gr-z  ~  O  a  Oq  ~  2  Qt  =  2*  +■  J ?izbfc  (x) 


6-=%e 


RESPONSE  NEAR  SOLID  SURFACE  TC  VARIABLE  A  M3 1  ENT  TEMPERATURE 
Analytical  Solution: 

The  following  problem  is  solved  in  ref.  2,  p.  297. 

’’The  region  x>0.  Initial  temperature  f(x).  Radiation 
(Newtonian  cooling)  at  the  surface  into  medium  at  ^(t)." 


T  =  h  ♦  I2  +  h 


where: 


_  r,  r  -tx-*'**  -(*+*') 

WihfL*  ™  -  * 


(57) 

(58) 

(59) 


MV 

snd:  I,=j[s£gL)  -ieA(t 

Each  of  the  above  three  integrals  will  now  be  expressed  in 


terms  of  the  functions  F„  and  G„. 

n  n 

Evaluation  of  I 


1* 

The  first  of  these  Integrals  is  handled  in  precisely  the  same 

manner  as  was  used  to  evaluate  T(0,t)  in  eqs.  3  to  9-  The  result 

is:  _ _  __ 

~ZV Jf? 

r  ivnt  rxi/ Lr^ytf/  -  z&VxtJ  J*  <*i\  , 


41 


cLx'  (61) 
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Evaluation  of  Igt 

The  second  integral  can  be  rewritten  as: 

Ii=-  AjF0(x+x',  t,A)  f(x‘)  dx' 

=  \f?(x)  iiF</dx')  d*' 

=  *'(*>  d*J 

0r!  -&  [Zi%) F, (x,  t,  t)  +  {Pfy  £,  dxj 

Evaluation  of  ly.  0 


(62) 

(63) 

(64) 


(65) 


The  third  Integral  can  also  be  expressed  simply  in  terms 
of  the  "F"  functions.  Thus: 


I,  =■  xi  JF,  (x,  t-rt  4)  j>(r)  dr 

Or,  with  the  use  of  eq.  51,  one  obtains: 

h  =£j$Hr)dT  =  -i  [*£!+&) dr 

T 

=  AjsF  <!>(r)  +  &  J F  4>  '(r) dr 

=  / '%(*,&)& -  ffr  tl(T)dr 

~&[F,  fa)+-L  f  t‘(o)J  +  k  ff  (fCridr 


The  general  result  is: 

DETERMINATION  OF  THE  TEMPERATURE  COEFFICIENTS 


<5 Mo) 

K* 


#• 


<f>  (r)  dr 


(66) 

(67) 

(68) 

(69) 

(70) 

(71) 


The  final  result  for  T(x,t)  is,  exclusive  of  error  terms: 

TM-k  -t±,Fc°)FJx,U)  +J 

+  (72) 

The  derivatives  appearing  in  eq.  72  can  be  expressed  in  terms 

of  finite  differences  .  If,  at  time  zero,  the  space  distribu¬ 
tion  of  temperature  can  be  expressed  by  a  second-degree  poly¬ 
nomial  in  "x" ,  and  the  ambient  temperature  as  a  linear  function 


"Numerical  Calculus,"  by  W.  2.  Milne,  Princeton  University 
Press,  Princeton,  N.  J.,  1949* 


of  "tw,  the  following  expressions  apply  to  the  various  derivatives 


f(o)~Te*  ,  (73) 

f‘(o)  =  (+T-  3Ta t  T,)U»)_  (74) 

f(o)  -  (To  -  z  T,  +  Ta  )  (&x )  2  (75) 

=•  ~Q( *o)  (76) 

-  ft  (At  -o)  -  TX+o) }  Cat)'  (77) 


These  expressions  are  used  in  this  paper,  although  the  result 
contained  in  eq.  72  applies  to  polynomials  of  arbitrarily-hlgh 
degree. 

When  the  finite-difference  expressions  73"77  are  sub¬ 
stituted  into  eq.  72,  the  coefficients  of  the  various  equally- 
spaced  temperatures  can  be  assembled.  For  the  case  where 
t  =  At  and  x  =  J(  Ax)  these  coefficients  are  given  in  eqs.  28-32 
of  the  text.  In  presenting  these  coefficients,  it  is  convenient 
to  use  the  dimensionless  sequence  of  functions  defined  by: 

F*  =  Fn/(Ax)n 


(78) 


